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Abstract 

We present a method to compute the density of states induced by a finite 
density of non magnetic impurities in a d-wave superconductor, in the uni- 
tary limit of very strong scattering centers. For frequencies very small as 
compared to the superconducting gap (w <C Ao) the additional density of 
states has the leading divergence Sp(u) ~ rij / ^|2oj| In 2 |w/Ao|J. This result 
is non perturbative. 
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As a consequence of the breakdown of Anderson's theorerrJll, when impurity scattering 
violates the symmetry of the superconducting state, the superconducting energy gap is de- 
pleted and impurities act as strong pair breakers. This is the case in s-wave superconductors 
with magnetic impuritiesili for which there is creation of bound states in the gap. Since they 
are not in the same symmetry representation, non magnetic impurities act as pair breakers 
in unconventional superconductors with higher orbital momentum such as d-wave super- 
conductors. As a result, the measured low temperature properties of YBaCu2O70 displays 
a remarkable sensitivity to the presence of non magnetic impurities: the critical temper- 
ature, for example, is suppressed even for very small density of impurities!. Furthermore 
the d-wave superconductor is special, due to the presence of gap nodes which prevents the 
complete freezing of scattering processes at low energy. 

The standard method to treat these problems of disorderly combines the T-matrix ap- 
proximation with standard impurity averaging techniques. For three dimensional systems 
such as polar superconductors or heavy Fermion superconductors^ the standard pertur- 
bative approach is reliable. In the limit of low impurity concentration rk. a perturbative 
expansion in n« leads to a finite density of states at the chemical potentialQO. 

For two dimensional systems (to which it is believed the high-Tc cuprates belong) the 
standard procedure of averaging over impurities may be complicated by the appearance of 
logarithmic singularities in the perturbative expansion of the single electron self-energy@. 
Such a situation appears in a variety of two dimensional systems characterized by a Dirac- 
like canonical spectrum. This remark has cast some doubt upon the validity of perturbative 
expansions in n« because at each loop-level logarithmic divergences prevent the series to 
converge. Nevertheless self-consistent versionsli^rEi of the standard averaging technique have 
been performed showing a finite density of states at zero energy. On the other hand, some 
non perturbative methods have been used to treat a weak disorder potentiafi'BHlI. The 
solution then depends on the symmetry of the pure system. In the special case of the d-wave 
superconductor the density of states still vanishes at low energy in the presence of disorder. 
Recently Ref. [L^ concluded that the density of states is finite above a very low energy scale 
(essentially the level spacingof a localization volume), below which a pseudo-gap appears. 
Early numerical simulationsEj seem to confirm finite density of states at zero energy. 

The issue of finite density of states is crucial for the conduction properties in the dis- 
ordered compound. Indeed, if there exist states in the gap, possible anomalous overlaps 
between well separated impurities can induce a new conduction mechanism entirely through 
impurity wave functions in a so-called 'impurity band^S and may lead to derealization. 
On the other hand, the results of Ref. [18] concluded that the states are localized. 

Here, we reexamine the issue of the density of states in a dirty d-wave superconductor. 
We consider the limiting case of a dilute concentration of identical impurities in the unitary 
limit. In other words, each impurity is a strong s-wave scatterer which is represented by 
an infinitely strong point-like repulsive potential (infinite scattering potential Vo) whose 
position is random. Using standard perturbative techniques, .previous studies concluded 
that the density of states was finite at the Fermi energyQoB'tallil. In contrast, using non- 
perturbative techniques, we find that the density of states is singular at the Fermi energy, 
p(u) = nil (2|c<j| In 2 |c<j/Ao|). The origin of this singular density of states is the existence of 
impurity bound states at the Fermi energy (E F = 0). It has been shownffl that one single 
impurity in the unitary limit creates a bound state at E — 0, with a spatial envelope that 



2 



decays as l/(R\nR). For many impurities these states overlap, but our result indicates that 
a singular density of states remains at zero energy of the form p(oS) = rii/ (2\u>\ In 2 |u;/A |). 

In this paper, we introduce a new method to calculate the leading divergence in the 
density of states induced by N non magnetic impurities in a system of two dimensional 
Dirac fermions in the unitary limit. We will outline the general features of the proof and 
present a different derivation than in our previous worked. The question of averaging will 
be discussed in details and we show how the result is exact provide the average over the 
different configurations of impurities is performed at the end of the calculation. Our result 
can be applied to several systems of Dirac fermions in two dimensions. We explicitly treat 
the d-wave superconductor. 

The paper is organized as follows. In section I we introduce the model of BCS d-wave 
superconductor, diagonalize it and define notations. In section II we establish the T-matrix 
equation starting from the equations of motion. We also introduce the matrix M as the 
inverse of the T-matrix. In section III we prove a sum-rule useful to calculate the density 
of states. Section IV is dedicated to a detailed study of the structure and matrix elements 
of M. Section V contains the heart of the proof. We first discuss the effect of taking the 
unitary limit. Then we indicate which are the quantities where the divergence appears, 
leading to the singular density of state that we found. In the conclusion we discuss how our 
result relates to the different theories of disordered superconductors. 

I. THE MODEL 

The generic Hamiltonian for a d-wave superconductor can be written 

#o = E4[ £ k <r 3 + AkaJ k . (1) 

k 

It describes BCS quasiparticles with the kinetic energy e k = W (cos k x + cos k y ) — /i (/i is 
the chemical potential) in the presence of the spin singlet superconducting order parameter 
Ak = Ao (cosfcz — cos k y ). Distances are measured in units of the lattice constant. The Oi 
are the Pauli matrices in the particle-hole space 




The spinor 0^ = \cl_p c_k,jj creates a particle and a hole with momenta k and — k, respec- 
tively. We shall present our results in terms of a d x 2_ y 2 state even though our conclusions 
apply more generally to any state where Ak vanishes linearly along a direction parallel to 
the Fermi surface. 

Instead of using the Nambu formalism, we work with the diagonalized version of ([!]) in 
order to access directly the properties of quasiparticles. The Bogoliubov transformation that 
diagonalizes H is given by 

Ck| = ^k«k - ^kAc , "k = ^kCkt + ^kcl k | > 
cL k | = ^k«k + «kAc , A< = -^kCkT + WkcL k |> (2) 
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where a k and /3 k create a particle and a hole with momentum k and the coefficients and 
t> k satisfy 



1/2 ( 1 + ^ 



Mk^k = 7; , 



with cj k = yej; + Aj;. Given the short-hand notation 



the BCS Hamiltonian can now be rewritten 



Ho = e E ^(-lriMk," • 

k f=0,l 

The disorder is introduced through N repulsive scalar potentials V located at random 
positions in the lattice: 

N 

Hi = V J2 E c tcia ■ (5) 

i=lff=U 

The full BCS Hamiltonian H = H + Hj describes a dirty d-wave superconductor. 



II. THE T-MATRIX EQUATION 
A. The Hamiltonian 

With the help of a Fourier transformation to the reciprocal lattice, the impurity potential 
becomes 

Hi = ^E E e*( k - k ')- R * (c kT c k , T + c^cvi) , (6) 

where V is the volume of the system. Rewriting the impurity term in terms of quasiparticles 
gives 

c kT c k'T = «kWk'a k «k' - ^kWk'A^k' - «k^k'a k A<' + v^v^ Ptf3k> ■ 

Thus 

i 

c icT c k'T = E (-iH-ir'^WV'LV'k'z/ , (7) 

i>,i/=0 

where we have introduced the short-hand notation 



tk,0 = 



Similarly, 



3 i(k-k')-Ri 



i k,k',j/,i/' 



^— *k^kV — t-kis+it-k'v'+l V'ki/V'k'i/' 



(8) 



cl k | C -k'| = - (^k^k'"k a k' + W k U k '/5 k a k' + Mk^k'CtjA' + Mk«k'/4 Ac') + COnst. , 

and, neglecting the constant term, we get 

c -ki c -k'i = -^WiW+iV'Lvw • (9) 

In summary, the random BCS Hamiltonian can be written 



(10) 



B. Equations of motion 

As the impurities break translation invariance, the anomalous two-point function 



GT q (T) 



^(r)V^(O) 



depends on two momenta. The equations of motion are 

Mcq^qk' - Okk'<W 



(11) 



(12) 



q,m 



where ££ k / i s 



££k' = [9r + (-l)^k]<W<W 

+ ^E^ (k - k,) - Ri (-i)^(-i)^k,w 



V 



Ei(k-k')-Riy. y. 



(13) 



We are dealing with a problem of non-interacting particles scattered by the static potential 
generated by N impurities. As we shall see, the anomalous Green function in Eq. (|TT| ) can 
be solved by inverting a 2N x 2N matrix. To this end, define the one-point functions 



Gl(r) 



d T + (-l)^k 



(14) 



(15) 
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together with 



A: 



B 



ij 



d 2 k 
d 2 k 

WT 2 



e - 4 k.(R,-R,) ( _ 1) n Ukv+iG i/ 



(16) 
(17) 
(18) 



The integration over k is to be performed over the first Brillouin zone. Equation (|i~3"D is 
inserted into the equations of motion ([121) . We then solve a 2N x 2N system of linear 
equations (see appendix |A|). This gives 



G kk> - G kk> ~ N -k^ ■ M 1 ■ N k v 
with a vector made of the 2N components 

/ glMi) \ 



(19) 



N 1 



N 2 



kv 



(20) 



and M is a 2iV x 2N matrix defined by 

M = 



I + V A V Q B 
V B I + V C 



(21) 



with A, B, C are N x N matrices whose matrix elements are respectively Aij, Bij and Cy. 
/ is the identity matrix. With these definitions the T-matrix equation can be written 



(22) 



with f = -QM- 1 . 



III. A SUM RULE FOR THE DENSITY OF STATES 



The increment in the density of states induced by the impurities is 

1 



5p(u) = Im J2 5G&(u} + i0 H 



7T 



(23) 



kv 



where 

6G& = -f ■ ■ N kv . (24) 
Recall that the summation over k is restricted to the first Brillouin zone. We can rewrite 
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'■j 



V 



If we go to Matsubara frequencies and define 



(25) 



^kv 



we notice that 



E (-VoNLNi 

ku 



ku 



-l)^ k ' 



d * 



and then 



J26G 



kk\ lUJ n 



ku 



—Tr 
V 



M 



.1 dM 
diui n 



(26) 



where the Trace is assumed to run over k, v as well as over the matrix indices. By convention 
we call 



5G = £ 5G& . 



(27) 



kv 



The equation ( 26]) can also be written (In Det = Tr In) 

1 dflnDetM 
5G(ia; n ) = — 



(28) 

We used this expression in a previous paperEI in order to derive the additional density of 
states. In this paper we will prefer to use the following expression 



5G{iu n ) = — Tr 



M~ 



dM 2 

2 diuj r , 



(29) 



IV. THE MATRIX ELEMENTS OF M 

As seen previously the evaluation of the density of states relies on calculating the deter- 
minant of M. We begin by evaluating its matrix elements. In this section, we just quote the 
result; the explicit calculation being given in Appendix ||. In order to perform this calcula- 
tion we make the assumption that the energy uj Ao. This will enable us to linearize the 
spectrum for small energies. Second we assume that W = Aq and the chemical potential 
fj, = in order to simplify the calculation. We will see later what happens when these two 
conditions are relaxed. Under these conditions there are four nodes in the Brillouin zone 
(cf. figure 0) located at (±|, ±|). With k' = (f , f ) + k we get 
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FIG. 1. Linearization of the energy spectrum at very low energy. In the first Brilloin zone we 
have four nodes centered at (±7r/2, ±7t/2) around which the energy spectrum is linearized. There 
are some symmetry relations between the four sud-zones called BZ 1, BZ 2, BZ 3 and BZ 4. The 
dotted line represents the points where £k = cos k x + cos k y = 0. The linearized spectrum around 
the nodes is represented in a third dimension in the plane (k x , k y ). 



LUI = W 2 



(cos k' x + cos k' y — fx) + (cos k' x — cos k'^j 



2W 2 k 2 . 



(30) 



Thus, co>k 
we find 



Dk with D = y/2W. after integrating over the four nodes in the Brillouin zone, 



where 



whereby 



and K Q is 
Appendix I 



Aij = A°(Hij) + ^(Ry) 
Bij = B 1 (Hij) , 



iuj n ) 2 - col 



ik R 



RiOr, 



D 



F (R) = 2 cos f (74 + Ry) + 2 cos §(i? x - i^) 

the Bessel function of rank zero. Note that \u n \ = 
1 Eq. 



(31) 



lUJ r 



(32) 
(33) 

) 2 . As shown in 
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-ik R 



k 

^(R) 



2 v / 2tt£) : 



with 



JF 1 (R) = 2 sin f (i^ + Ry) (cos + siny?) + 2 cos f (i^ — Ry)(cos<p — sin <p) 



(34) 



(35) 



where <p is the angle between R and the x-axis and K\ is the Bessel function of rank one. 
In the same manner 



ik R 



^ 2 (R) 



(Or. 



2 v / 2vr J D 2 



Ruo r 



D 



(36) 



with 



JF 2 (R) = 2 sin %(R X + R y )(cos <p — sin tp) + 2 cos f (i£c — R y )(cosip + sin</?) . (37) 



Note that the point R = is rather special with A°(0) = f^gf In |w n /£>| and with ^(0) 



V. EVALUATION OF THE DENSITY OF STATES 



A. Unitary limit and low energies 

In this section we will evaluate the leading term in the density of states in the limit of 
low frequencies. When < i?|a> n | <C D we have 



K 



R\uj n \ 



In 



R\io n \ 



D J \ D 
so that in the limit \u n \/D — > we have 

A°(R) 



R\u> n \ 

D 



D 



\w n \R 



(38) 



(39) 



In the limit of low frequencies and for 0, this enables us to neglect ^4°(R) as compared 
to A 1 (R) in the evaluation of the matrix M in Eq. ||. For R = 0, ^t°(0) and A l (0) are 
negligible as compared to ^(R) for R ^ 0. In the sequel we can safely avoid the point 
R = in summations over R that occur during the evaluation of M 2 . Second we notice 
that the Bessel function Kq and K\ have an exponential cut-off at R ma x = D/\u n \ so that 
we can safely use the approximation: 



A 1 (R ij ) ~ { 2V2ttDR 
, 



Kl 



if R < D/\u r , 
elsewhere , 



(40) 
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^ (R " ) if R<D/\u„ 



B 1 (R lj )~\ 2^D Rij ' (41) 
, elsewhere . 

An important remark to make is that the a;-dependence of the matrix elements of M appears 
only through the upper cut-off of the Bessel functions. In what follows, we make the crucial 
assumption of the unitary limit, i.e., that Vq — > oo. Recalling the form of M in equation fl2~T|) 
we see that in this limit the identity matrix in M becomes negligible when compared to 
and Cij. 



B. The divergences appear 



Our aim is now to factorize the leading divergences in this problem. In order to make 
the divergence apparent, it is more convenient to work with 



5G(iu n ) = Tr 



M~ 



dM 2 



2 diu) ri 



(42) 



From section |V A| , we expect logarithmic factors \n\D/uj n \ to appear. An important point 
to stress is that for any configuration of the impurities, we will always find some factors 
In \ D/u n \ in M 2 . Within the unitary approximation we have 



M 



A 2 + B 2 AB-BA 
AB-BA A 2 + B 2 



(43) 



To see that In \D/u n \ is necessarily present in the diagonal terms of M 2 , we estimate 

1 [{J*f + {F 2 ) 2 } 



E 



(2tt£)) : 



p2 



(44) 



where the summation over j is restricted to < Rij < D/\u n \. Provided the impurities are 
rather homogeneously scattered in the system (around each impurity site Rj one can find 
a macroscopic amount of impurities inside a circle of radius D/\u n \, cf. figure || ), we can 
take the continuous limit, 



(45) 



Now [T ) + [T ) are oscillatory but always positive so for all impurity at site Rj we have 



Mf i = C\n 



D 



(46) 



where C is a constant. 
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FIG. 2. The point Rj and Rj where the impurities are located are represented in this figure. 
The integrand in Eq. ( f45| ) centered around the point Rj is nonzero only when Rij < D/\u n \. 



The biggest coefficients in M 2 are situated on the diagonal. To see it, we distinguish 
between the off-diagonal elements that are diagonal in the particle-hole grading and those 
that are not. The magnitude of off-diagonal elements that are diagonal in the particle- hole 
grading can be written 



Ml 



E 



• r ij' r jk 



• r ij- r jk 



(47) 



As i 7^ k, M 2 k picks up an oscillatory prefactor (a combination of e ±t7T ^ 2Rik e ±%ip as in 
eqns. (|33|), (|35|), (137|) ). Furthermore, the logarithmic divergence gets rescaled by Rik'- 



Ml 



DR 



ik 



C Osc(R ik ) In 
where |Osc(Rjjt)| < 1. Hence for all sites k 7^ i we have 

< 1 . 



(48) 



Ml 



M 2 



(49) 



In turn, the terms in the off-diagonal blocks of M 2 with respect to the particle-hole 
vanish. Indeed, if we denote the off-diagonal elements of M 2 with respect to the particle- 
hole grading by 



(50) 



we have 
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FIG. 3. Vanishing of the element Tj&. The summation zone is the surface of intersection of the 
two disks. The first term in ( |5l|) is represented in the upper drawing whereas the second one is 
represented in the lower drawing. For each summation point Rj in the upper drawing there is a 
symmetric one Rj in the lower drawing such that Ry = Rj^ and Ry = R,-fc. 



ik Y ( 27rD ) 2 \ ^j R 3k RijRjk ) ' 1 J 

where here the sum runs over the points R., such that Rij < D/\u n \ and Rj k < \u n \. Noticing 
the symmetry of JF 1 and T 2 under the transformation (p — > if + tt, we can show that the 
two terms on the r. h. s. of Eq. (^TJ) cancel identically. Indeed the integrands — j^ 13 - and 

JC "^' fc ' ) are represented respectively within each circle of figure |] ( both of these terms have 
a cut-off at Rmax = D/\u n \). The summation zone is the surface of intersection of the two 
disks. The first term in (|5l|) is represented in the upper drawing whereas the second one is 
represented in the lower drawing. For each summation point Rj in the upper drawing there 
is a symmetric one Rj in the lower drawing such that Ry = R^ and Ry = Rjfc. Thus 
T ik = 0. 

In conclusion we find it convenient to define a matrix S by 



M 2 = Cm 



D 



LO r , 



S , (52) 



where C is the constant defined in P6|), independent of disorder. The matrix S depends on 
the particular configuration of the impurities in the system. It satisfies 

\Sij\ < 1 • (53) 
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C. Asymptotic value of the density of states 

Substituting the value of M 2 from ( j5^) into the equation (f42|) we find 

8G(iu n ) = T div (iu n ) + TZ(iuj n ), 
where (■) denotes the average over disorder, / is the 2N x 2N density matrix, and 
1 d 



2V diu, 



lnln(£>/KI) Tr I 



\calR(iuj r: 



1 (TrS-'J-S) 



2V 



diu r 



(54) 



(55) 



Then the first term in Eq. ([54]) %a v is responsible for the singular density of states that we 
obtain. Indeed it gives rise to 



N 



(56) 



V iuj n In \u n /D\ ' 

where we recall that N is the number of impurities. After analytic continuation (remem- 



ber that \u n \ = 
negligible, we get 



[lUJr 



) 2 , Eq. (|B16j ) ) and assuming that the reminder 1Z in Eq. (|54]) is 



and thus 



5p(uj) 



■—riilm 



11; 



+ i8) ln(^) 



(57) 



1 



\U)\ 



ln 2 (| W |/ J D) + ( 7 r/2)^ 



(58) 



where rii is the density of impurities in the system. We note that this expression is normal- 
izable: 



D 



D 



5p(u))du) = 2rii . 



(59) 



In order to prove the result ( |58"D we still have to show that the reminder 1Z in fl4*2| ) is negligible 
as compared to Tdm- I n order to do this we have to give some insight about the form of S~ x . 



D. The form of S" 1 

The matrix S is invertible (since M is invertible) and we will find a reasonable candidate 
to the inverse of S in order to give an estimation of the reminder TZ. Inverting S means we 
can find a matrix S^ 1 such that for any given pair of sites i and k 

E Mr* 1 = 6 « ■ ( 6 °) 

3 

We introduce a pictorial representation of S by drawing a disk ( called S-disk) of radius 
\D/ui n \ centered at R.;. To each location Rj of an impurity there corresponds a matrix 
element Sij which depends on the vector R^- = Rj — Hj. From Eqs. (^0|,|4lD we recall that 
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(a) (b) 

FIG. 4. Volume of summation of J2j SijSj^ in two configurations. In (a) we have i = k and in 
(6) we have i ^ k but i still very close to k. 

\Sij\ < 1 , for Rij < \D/u n \ , . . 

Sii = , for i?y > IDKI , l ° j 

inside the disk, 5 has some non vanishing matrix elements \S{j\ < 1. Outside the disk, 
Sij = 0. It's important to note at this point that the only dependence on u in the S matrix 
comes from the cut-off. In order to satisfy fl60j) we must presume that S -1 has the same 
cut-off \D/u n \ as S. Hence we represent again by a disk (called S _1 -disk), but centered 
around this time. Now the summation J2j ^ijSj'j} runs over the intersection of the S-disk 

and the S _1 -disk. The key difficulty in order to invert S is to find a matrix S" 1 , such that 
when the two disks have the same center we have 

3 

whereas when the two centers differ, even by a small amount, we have 

J2j Sij&jk = ; k i . 

This is illustrated on figure |] where on the left side (case (a)) the two circles are centered 
at the same point and on the right side (case (b)) the two centers differ by a tiny amount. 
In both cases the intersecting area of the two disks is almost identical, but in case (a) the 
result has to be 1 whereas in case (b) it has to be 0. 

The matrix elements of S inside the disk are random, but the condition \Sij\ < 1 is 
independent of the realization of disorder. The worst possible situation for differentiating 
between cases (a) and (b) is when all the matrix elements inside the 5-disk have their 
maximum value 1. 

We define the external boundary of the S'" 1 -disk. By external boundary we mean the 
circle exactly adjacent to the disk and external to it. Upon the external boundary the 
matrix elements of S" _1 are defined as non vanishing and negative, so that the summation 
over it compensates the summation over the intersection of the S'-disk and the S' _1 -disk. In 
case (a) the external boundary doesn't touch the S'-disk and thus has no effect upon the 
summation over the intersection of the S'-disk and the S~ 1 -disk. Alternatively in case (b) 
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the summations over the intersection of the 5-disk and the S _1 -disk and over the external 
boundary cancel out. 

Let's take an explicit example where = 1 if Rij < \D/u n \ and Sij = elsewhere. 
We define S', 1 in the following way. For < \D/u n \, S^ 1 is proportional to a random 



configuration of ±1 such that after integration over a disk of volume V = n \ D/u n \ we get 

V 

D 



(62) 



By the central limit theorem, there are many random configurations of ±1 that verify this 
condition. We then take the external boundary to be proportional to with the 

same proportionality constant. Thus 



1 

sr, 1 = -a 



if < \D/u n \ 



if R, 



7T 



\D/u n \ 



(63) 



, if R i:j > \D/uj n \ , 

where A is a constant. We notice that when and R& are infinitely close, then 



boundary 



D 



and exactly compensates the summation over the volume inside. The proportionality con- 
stant A is fixed so that 



In summary 



thus A 



UJr, 







q-l 




nD 


A-l 


= , 



;±1) , if Rij < \D/w ri 



if Ri 



\D/oo n \ 



(64) 



if R^ > \D/u n \ . 



Now for each intermediate case where k ^ i (cf. figure |5]) we want to be sure that the 
intersection of the S*-disk and the S _1 -disk compensates the sum over the external boundary 
of S 1-1 which crosses the S'-disk. This is obviously not the case for any configuration of 
random ±1 in S' 1 but we believe there is one (and actually only one because there is only one 
inverse for SI) configuration which satisfies it for all positions of Rj and R&. As represented 
in figure [VD| we call 9 the angle made by the center i with the two points A and B, 
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FIG. 5. Compensation between the volume and boundary. 



intersection between the two circles. The area of summation is given by A e = 9R 2 — AR 2 sin 9 
and the intersecting arc's length by L e = 9R. They both scale in respectively R 2 and R 
with varying prefactors. So the desired configuration of random ±1 in the S'~ 1 -disk has to 
be "denser" towards the center of the circle than towards the boundary. We are not able to 
write down explicitly this configuration of random ±1 inside the area of the S^-disk, but 
actually it doesn't matter, because as we will see in the next paragraph the evaluation of 
the reminder doesn't depend on it. 




E. Evaluation of the reminder TZ 

Now we evaluate 

An important point is that as S depends on uj n only via its boundary, dS /duj n is a matrix 
with non zero values only on the external boundary of the S-disk. explicitly, taking our 
example where Sij = 1 if Rij < D/\uj n \ we get 
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dS 

duj n 

dS 







if Rij = D/\u n \ , 



elsewhere . 



(66) 



Since the matrix elements in the S 1 -disk have random (positive and negative) signs whereas 
upon the external boundary they have constant sign (negative here), the maximum value 



of 



S^dS/dUr 



is reached when i = k, that is when the external boundary of the S'-disk 

(where the matrix elements of dS/du n are non vanishing) matches the external boundary 
of the S^-disk (cf. Fig. g). 
In our special case we get 



S 



_! dS 



du) r , 




(67) 



Since 
equal 



dS/du r , 



has its maximum value when all the matrix elements of S inside S'-disk 
we have indeed evaluated an upper bound for the reminder TZ. The reminder 1Z 
is thus negligible compared to the leading divergence in the density of states in the limit 
\u n \/D < 1. 



VI. DISCUSSION 

The first question to address is what happens when the different conditions under which 
our calculation was performed are relaxed. First consider the more realistic situation where 
the bandwidth W and the superconducting gap A are not equal (experimentally, we have 
A — 0.10 W). According to Ref. ^5] the power law dependence of the matrix elements of M 
(cf. equations fl^O] ) and (0)) is still preserved, but gets an overall prefactor of A /W. The 
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scale under which our calculation is valid is now \u\ <C VoA /W. In addition, the hopping 
matrix M will show a strong spatial anisotropy0. This anisotropy can be absorbed into the 
overall prefactors JF 1 and T 2 (resp. eqns. (|33|), fl35|) and fl3T|)) entering the definition 
of the matrix elements Ay, Cy and S^. Since only the square of these factors enter the 
leading divergence (cf. Eq. flH|)), we believe the result stays unchanged. 

What happens if the bandwidth are not symmetric anymore, that is if ji 7^ but still 
/i< A ? 

First the nodes are moved away from the point (±7r/2, ±7r/2) so that transversal nodes 
are now separated by the vectors Q = — 5), 7r(l — 5)) and Q* = (— 7r(l — 5), 7r(l — 5)) 
where 5 = /i/A and /z is the increase in the chemical potential. This leads to a change of 
the phase factors in Aij and Namely we get JF°(i?) = 2 [cos (Q • R/2) + cos (Q* • R/2)]; 

T 1 (R) = 2 [sin (Q • R/2) (cos (p + sin (p) , 

— sin (Q* • R/2) (cos ip — sin <p)] , 

and 

T 2 {R) = 2 [sin (Q • R/2) (cos up - sin ip) , 

— sin (Q* • R/2) (cos <p + sin <p)] . 

As 5 is a small parameter, this change in the phase won't affect the existence of the logarith- 
mic divergence in M 2 . Additionally, away from half filling, the bands of quasi particles and 
quasi holes become asymmetric to account for the removing of particles in the system. The 
difference induced in A^, Bij and CV,- comes from the highest part of the energy spectrum, 
where k ~ 1 and uo ~ D and shouldn't affect our result. 

Our solution is valid under the assumption of unitary limit, meaning that Vq is the largest 
scale in the problem. As shown i to non interacting single impurities is associated 

the creation of a bound states decaying as l/(i?lni?). Following Balatsky et al. (@) this 
would lead to derealization due to the formation of an impurity band. The result we find for 
the density of states is indeed reminiscent of a Dyson-like singularity ( in d=l, a Dyson-like 
singularity corresponds to a density of state diverging in 1/ (Jc<j| ln 3 (|c<j|/.D)) ), associated in 
one dimension with delocalizationBZl. 

What happens when the condition of unitarity is relaxed is still very much an open 
problem. In the case of weak disorder, some a-model analysis have been performedBtllll 
concluding to a vanishing density of state under an energy scale E2 = Dp/C, 2 , where Dp is 
the bare diffusion constant and £ is the localization length. This result is strongly supported 
by symmetry considerations. Indeed, a disordered d-wave superconductor belongs to class 
Cj, according to the classification of ref.il, meaning that the Hamiltonian is invariant under 
time-reversal symmetry as well as spin rotation symmetry. According to random matrix 
theory a universal behavior is expected, inducing a vanishing density of states on the scale 
of the level spacing induced by the finite localazation length. 

Consider one impurity with scattering potential Vq. The effective potential at the im- 
purity site can be evaluated exactlyili! and is given by V = i_v \cj\in\D/u\ • ^ n ^ ne un itary 
limit ( Vq — > 00) we get V = , , ^i d / u i ■ This effective potential diverges when uj goes to zero. 
On the other hand, we notice that the derivation of non linear sigma models for disordered 
systems requires Gaussian disorder, and especially require that the average effective disorder 
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potential vanishes (V) = (the second moment is nonzero). If (V) is nonzero but is a con- 
stant as a function of energy, it can be absorbed as a redefinition of the chemical potential, 
but the case where the effective potential would diverge as u> goes to zero belongs to another 
universality class: the energy of an effective non-linear sigma model would renormalized to 
uj — S(co>), where S(a>) diverges when omega goes to zero. The energy scale under which the 
non-linear sigma model describes the diffusive modes then becomes inaccessible since the 
effective energy never gets close to zero. In our problem the result obtained on the density 
of states indicates that the self-energy is of the form H(u) = , - ] — r , diverging as u 

goes to zero, but still different from the one impurity case. The feynmann diagrams leading 
to such this self-energy will be studied in a future publication!!. We believe the method 
presented here, using the T-matrix equation takes care in a non perturbative way of the 
leading divergence in the unitary limit. One possible scenario which would reconciliate the 
two limits of weak and strong disorder is that the unitarity limit fixed point is unstable (as 
soon as Vq becomes finite, the effective potential saturates), but the cross-over regime close 
to unitarity is still very much influenced by the strong disorder fixed point. 
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the resonant density of states alone was not sufficient to induce delocalisation. 
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APPENDIX A: DERIVATION OF THE T-MATRIX EQUATION 



Starting form equation fll2"D , replacing the lagrangian of (13) into it and multiplicating 
from the left by G° gives 



i 

+ E e * '^Wfl^ki^kV = G^iAk'^' 



(Al) 



where 



Gl /"D \ _ v-v — tq-Ri / i\mj fimv 

kl/V"V — l^qm e I X J ^qm^qk ' 

LT k!/l rv V — Z^qm e V+1 U ' 



mv 
qk • 



(A2) 



Now we have two unknown functions G 1 and G 2 that we evaluate using formula (|A1~1): 
- + § E E e iq - (R <- R ^V) 2 G qm GUR*) 

+ ^EE^ (Rl+Rj) (-l) m W + iG' qm GL(R,) = (-lTe- ik - R H kv Gl ; (A3) 



i qm 



GURj) + ^EE eiq ' (Ri " Rj) (-l) m iqmV + lG' q mGL(R i ) 
" i qm 



i qm 



These two equations can be rewritten matricially as 

{-5 t3 + V A tJ ) Gl(R,) + VoB^GURj) = A^R,) , 
(5jj + V Cij) G ku (Rj) + VoBijG ku (Rj) = N ku (Rj) . 



Define the 2N vector V 



V = 



V 1 
V 



/ GL(Ri) \ 



2 



V 1 



ku 



V 2 



k;/ 



and the equation (|A5|) can be written 

MV k , = N k 



/ GL(Ri) \ 



(A4) 



(A5) 



(A6) 



(A7) 



But then equation (|A1|) becomes 



^kk' + ^jjN kl/ • V ki , — G^vSkk'tii'v' 



Insertion of (|A7Q yields 



G kk> - G kk> - N -k^ ■ M 1 ■ N k v 



(A8) 



(A9) 
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APPENDIX B: CALCULATION OF THE MATRIX ELEMENTS A u , B u AND Cjj 



„ ■ (2tt) 

d 2 k 



2 

ik-R, 



ul 



+ 



4 



(2tt) 2 \iu n — iuj n J rUJ] i/ 

where the integration runs over the first Brillouin zone. Since 



we get 



with 



Similarly 



and 



u 



k — 2 
1 

2 



£(1 + ^ 



-4 (Rij) — (2^)2 ; ,J_ ,J ■ 



k 



A 1 (R ij ) =f (f k ,i ( . % a e~» ' 



B 



d 2 k A 



" ' (2vr)2 (^ n )2 _ ^ 



k g-ik-Rjj 



(Bl) 
(B2) 



(B3) 



(B4) 



(B5) 



(B6) 



(B7) 



1. Evaluation of.4 (R,i) 
As we can see on figure [l] the spectrum has four nodes at the points 

P /7T7T\ JD { IT TT\ p I TT 7T\ p 

■Tl — {01 0)1 ^ 2 — \— 9, 9 ) , -T 3 — \— 9 , — 9 J ) -T 4 



7T 7T^ 

2' 2' 



:bs) 



Under the assumptions fi = and A = W, we can linearize the spectrum around each node 
in the following way: 



k' = (f,f)+k 
and we get 



W 2 



'cos k' x + cos k'y) + (cos k' x — cos k y 



I \2 



I ~D 2 k 2 , with D = V2W 



(B9) 



(BIO) 
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Similarly, 

e' k ~ W(k x + k y ) , ^ B1 ^ 



~ (k x - ky)A . 

In order to evaluate ^4°(R) we divide the integral into a sum of four integrals around 
each node : 

e -ik'-R, e -ik'R 

e -ik'R e -ik'R 

+ E s 2 2 + E """ T-TTTa ,,2 ' ( B12 ) 

In each term, the fact of developing around a particular node gives a specific prefactor, so 
that we have 

with 

^O(R) = e -»f(-R*+^) + e -*f (H»--Rv) _|_ e *f _|_ e i^(R x -R y ) _ (B14) 

Thus 

^°(R) = 2 cos + 2 cos -(i^-i^) . (B15) 

Now calling 9 the angle between k and R (cf. figure [7]), 

1 k dk e -ikRcos8 

d9- 



o (27r) 2 y (iu n ) 2 -D 2 k 2 



--^ U (R)^^7 / k dk 



2nD 2 Jo (oo n /D) 2 + k 2 



-T\B)-^-K {R\uj n \/D) , (B16) 



1uJ 



2nD 

where K is the Bessel function of rank zero. Note that we have defined \uj n \ = J —{iuj 



2. Calculation of ^(R) 
As previously, we can decompose the summation in the Brillouin zone into four parts: 



_ ,,-ik'R _ ,,-ik'R 

+ Y. % 2 + E 7^ 5 i ■ PIT) 
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FIG. 7. Angles between R and x and between k and R. 



If we call 6 the angle between k and R and tp the angle between R and the x-axis as 
represented on figure [7|, the first term in ( |B17| ) can be written 



e k > e 



ik'R 



k'eBZ! 



[ILOr, 



OJt, 



2^ 



d0 / fc dk 



D (k x + k y )e 



-ikR cos 







v/2(2tt) 2 (iw») 2 " -D 2 ^ 2 



(B18) 



Using 



/c x = fc(cos cos 9? — sin 6 1 sin <p) , 
fcj, = /c(cos 9 sin + sin sin </?) , 



(B19) 



we get 



Ti = e'^^+^cosy^ + siny?) 



v/2(2tt 2 ) V 

iD r 1 



1 7 2 j, /" 27r .„ ( C0S 61 + Sin 

k dk / 



-ikR cos 8 



e 2 



if (^+^)( COS(j0 _|_ s i nv9 ) _L — /" k 2 dk 

\ Tin h 



J x {kR) 



( /, ; '„) 2 - D 2 & 2 

—ikR cos 



2V2vr 1 



v/22 



D 2 k 2 



(B20) 



Thus, after summing over the four nodes we get, 



with 



7T 7T 

sin — (i^a; + By) (cos y? + sin if) + sin — (i?^. — R y ) (cos y? — sin if) 



JF 1 = 2 

The evaluation of i3 x (R) is done in the same way as the one of ^4 1 (R). 



(B21) 



(B22) 
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